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ABSTRACT 

We present new models of the evolution and dissolution of star clusters evolving under the 
combined influence of internal relaxation and external tidal fields, using the anisotropic 
gaseous model based on the Fokker-Planck approximation, and a new escaper loss cone 
model. This model borrows ideas from loss cones of stellar distributions near massive black 
holes, and describes physical processes related to escaping stars by a simple model based on 
two timescales and a diffusion process. We compare our results with those of direct TV -body 
models and of direct numerical solutions of the orbit-averaged Fokker-Planck equation. For 
this comparative study we limit ourselves to idealized single point mass star clusters, in order 
to present a detailed study of the physical processes determining the rate of mass loss, core 
collapse and other features of the system's evolution. With the positive results of our study the 
path is now open in the future to use the computationally efficient gaseous models for future 
studies with more realism (mass spectrum, stellar evolution). 

Key words: gravitation - methods: numerical - celestial mechanics, stellar dynamics - glob- 
ular clusters: general 



1 INTRODUCTION 

Dynamical modelling of globular clusters and other collisional stel- 
lar systems (like galactic nuclei, rich open clusters, and rich galaxy 
clusters) still poses a considerable challenge for both theory and 
computational requirements (in hardware and software). On the 
theoretical side the validity of certain assumptions used in statis- 
tical modelling based on the Fokker-Planck (henceforth FP) and 
other approximations is poorly known. Stochastic noise in a dis- 
crete TV-body system and the impossibility to directly model real- 
istic particle numbers with the presently available hardware, are a 
considerable challenge for the computational side. 

Detailed comparisons of the results obtained with the different 
methods for single mass isolated star clusters have been performed 
(Giersz & Heggie 1994a,b, Giersz & Spurzem 1994, Spurzem & 
Aarseth 1996, Spurzem 1996). They include theoretical models 
such as the direct numerical solution of the orbit-averaged ID 
FP equation for isotropic systems (Cohn 1980), isotropic (Heggie 
1984) and anisotropic gaseous models (henceforth AGM) (Louis 
& Spurzem 1991, Spurzem 1994) and direct TV-body simulations 
using standard iV-body codes (NBODY5, Aarseth 1985, Spurzem 
& Aarseth 1996; NBODY2, Makino & Aarseth 1992; NBODY4, 
Makino 1996; NBODY6++, Spurzem 1999, Aarseth 1999a,b, 
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2003) . All the cited work, however, only dealt with idealized single- 
mass models. There are very few attempts yet to extend the quan- 
titative comparisons to more realistic star clusters containing dif- 
ferent mass bins or even a continuous mass spectrum (Spurzem & 
Takahashi 1995, Giersz & Heggie 1996, Giirkan, Freitag & Rasio 

2004) . 

On the side of the FP models there have been two major recent 
developments. Takahashi (1995, 1996, 1997) has published new FP 
models for spherically symmetric star clusters, based on the numer- 
ical solution of the orbit-averaged 2D FP equation (solving the FP 
equation for the distribution / = f(E, J 2 ) as a function of en- 
ergy and angular momentum, on an (E, J 2 )-mesh). Drukier et al. 
(1999) have published results from another 2D FP code based on 
the original Cohn (1979) code. In such 2D FP models anisotropy, 
i.e. the possible difference between radial and tangential velocity 
dispersions in spherical clusters is taken into account. Although 
the late, self-similar stages of core collapse are not affected very 
much by anisotropy (Louis & Spurzem 1991), intermediate and 
outer zones of globular clusters, say outside roughly the Lagrangian 
radius containing 30 % of the total mass, do exhibit fair amounts 
of anisotropy, in theoretical model simulations as well as accord- 
ing to parameterized model fits (Lupton, Gunn & Griffin 1987). In 
contrast to AGMs the 2D FP models contain less inherent model 
approximations; they do not assume a certain form of the heat con- 
ductivity and closure relations between the third order moments as 
in the case of AGM. Furthermore, the latter contains a numerical 
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constant A (Spurzem 1996), which is of order unity, but its numer- 
ical value has to be determined from comparisons with proper FP 
or A-body models. 

Secondly, another 2D FP model has been worked out re- 
cently for the case of axisymmetric rotating star clusters (Einsel 
& Spurzem 1999, Kim et al. 2002, Kim, Lee & Spurzem 2004, 
Fiestas, Spurzem & Kim 2005). Here, the distribution function is 
assumed to be a function of energy E and the z-component of an- 
gular momentum J z only; a possible dependence of the distribution 
function on a third integral is neglected. As in the spherically sym- 
metric case the neglecting of an integral of motion is equivalent to 
the assumption of isotropy, here between the velocity dispersions in 
the meridional plane (r and z directions); anisotropy between that 
velocity dispersion and that in the equatorial plane ((^-direction), 
however, is included. 

Thirdly, there is an elegant alternative way to generate mod- 
els of star clusters, which can correctly reproduce the stochastic 
features of real star clusters, but without really integrating all or- 
bits directly as in an A-body simulation. They rely on the FP ap- 
proximation and (hitherto) spherical symmetry, but their data struc- 
ture is very similar to an A-body model. These so-called Monte 
Carlo models were recently redeveloped by Giersz (1996, 1998, 
2001), and by Rasio and collaborators (Joshi, Rasio & Portegies 
Zwart 2000, Watters, Joshi & Rasio 2000, Joshi, Nave & Rasio 
2001, Fregeau et al. 2003, Giirkan, Freitag & Rasio 2004). For an- 
other approach, reviving Henon's superstar method, compare the 
work by Freitag (Freitag 2000, Freitag & Benz 2001, 2002). The 
basic idea is to have pseudo-particles, whose orbital parameters 
are given in a smooth, self-consistent potential. However, their or- 
bital motion is not explicitly followed; to model interactions with 
other particles like two-body relaxation by distant encounters or 
strong interactions between binaries and field stars, a position of 
the particle in its orbit and further free parameters of the individ- 
ual encounter are picked from an appropriate distribution by using 
random numbers. A hybrid variant of the Monte Carlo technique 
combined with a gaseous model has been proposed by Spurzem 
& Giersz (1996), and applied to systems with a large number of 
primordial binaries by Giersz & Spurzem (2000) and Giersz & 
Spurzem (2003). The hybrid method uses a Monte Carlo model 
for binaries or any other object for which a statistical description, 
as used by the gaseous model, is not appropriate, due to small num- 
bers of objects or unknown analytic cross sections for interaction 
processes. The method is particularly useful for investigating evo- 
lution of large stellar systems with realistic fraction of primordial 
binaries, but could also be used in future to include the build-up of 
massive stars and blue stragglers by stellar collisions, for example. 

In the present and near future a wealth of detailed data on 
globular clusters will become available by e.g. the Hubble Space 
Telescope and the new 8m-class terrestrial telescopes such as Gem- 
ini and the Very Large Telescope (VLT), for extragalactic as well 
as Milky Way clusters. These data cover luminosity functions and 
derived mass functions, color-magnitude diagrams, population and 
kinematical analysis, including binaries and compact stellar evo- 
lution remnants, detailed two-dimensional proper motion and ra- 
dial velocity data, and tidal tails spanning over arcs several degrees 
wide (Koch et al. 2004). With detailed observational data such as 
from King at al. (1998), Piotto & Zocalli (1999), Rubenstein & Bai- 
lyn (1999), Ibata et al. (1999), Piotto et al. (1999), Grillmair et al. 
(1999), Shara et al. (1998), Odenkirchen et al. (2001), Hansen et 
al. (2002), Richer et al. (2002) to mention only the few recently ap- 
peared papers), an easily reproducible reliable modelling becomes 
more important than before. For that purpose a few more ingredi- 



ents are urgently required in the models in addition to anisotropy 
and rotation: a mass spectrum, a tidal field, and the influence of 
stellar evolution on the dynamical evolution of the cluster. 

While it is easy in principle to include all these in a direct A- 
body simulation, and considerable effort goes in the construction of 
new hardware and software for that purpose (Hut & Makino 1999, 
Makino et al. 1997, Makino & Taiji 1998, Aarseth 1999a,b, 2003), 
the life span of globular clusters extends over tens to hundreds of 
thousands of crossing times, taken at the half-mass radius, which 
requires a high accuracy direct A-body code for its modelling. De- 
spite the enormous advances in hardware and software efficiency, 
still modelling, say, of a few hundred thousand particles is very 
tedious. For parameter studies one has to rely on lower particle 
numbers and prescriptions to scale the results to larger A (Aarseth 
& Heggie 1998, Baumgardt 2001). So we still need the fast but 
approximate theoretical models. Takahashi, Lee & Inagaki (1997) 
published first results for anisotropic clusters in a tidal field, and 
Takahashi & Lee (2000) extended their study to multimass clus- 
ters. Takahashi & Portegies Zwart (1998, 2000), Portegies Zwart 
& Takahashi (1999) compared the influence of stellar evolution- 
ary mass loss as measured in A-body and FP models. However, 
still the tidal boundary poses an unsolved problem in the pure point 
mass case. Another difficulty is that A-body models and theory are 
difficult to match, since an energy cutoff as it is usual in the FP 
models has not yet been applied in A-body simulations (with the 
exception of one of the models kindly provided by E. Kim for our 
comparisons, see below); usually direct A-body models employ a 
tidal radius cutoff. 

Moreover, in the theoretical models one has to decide which 
criteria to use for a star to escape. In contrast to isolated clusters 
energy and angular momentum are not exact integrals of motion in 
a tidal field, and one should consider more appropriate quantities 
like the Jacobian. However, even then it is possible that a star sat- 
isfying an escape criterion stays for many orbital times close to the 
cluster and even can be scattered back to the cluster (Baumgardt 
2001), so from an observers viewpoint they do not escape imme- 
diately; indeed Fukushige & Heggie (2000) propose a new form of 
an energy dependent escape time scale, which is formally infinite 
for stars at the tidal energy (note a similar problem in relation to 
dwarf spheroidal galaxies in tidal fields raised by Kroupa 1998 and 
Klessen & Kroupa 1998). 

In this situation we propose to look at simple cases and sim- 
plified models in more detail first. For example it appears to be 
useful to separate the dynamical effects of stellar evolution from 
those induced by the interplay of internal relaxation and a tidal 
field, and to look at a single mass case first. A gaseous model of 
a star cluster (here used an anisotropic gaseous model or AGM, see 
http://www.ari.uni-heidelberg.de/gaseous-modeli is a very sim- 
ple tool with which to understand, on the basis of an idealized 
model, the physical processes acting on star clusters. It has been 
very successfully used to detect gravothermal oscillations (Bet- 
twieser & Sugimoto 1984) and to discuss effects present in multi- 
mass star clusters (Spurzem & Takahashi 1995) and in systems with 
primordial binaries (Heggie & Aarseth 1992). 

Thus in this paper we present a new approach to include a 
tidal field into AGM of star clusters, which has not been tackled 
before. The processes of escape and relaxation are treated in a sim- 
plified parameterized way and compared to models using an orbit 
averaged Fokker-Planck equation and direct A-body simulations. 
Also the different parameters of our escaper model are varied to 
understand their influence on the results. We will demonstrate that 
the gaseous model helps in the physical understanding of the tidal 
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escape process and allows a deeper insight into the processes going 
on in other models. With this the way is paved for an application of 
this same gaseous model including more realistic properties of star 
clusters, such as a mass spectrum, stellar evolution and primordial 
binaries. 



2 THE MODELS 

The anisotropic gaseous model (AGM) is based on local moments 
of the stellar velocity distribution function. From the local Fokker- 
Planck equation (Rosenbluth, McDonald & Judd 1957) moment 
equations up to second order are taken and closed in third order 
by a phenomenological heat flux equation (Lynden-Bell & Eggle- 
ton 1980, Bettwieser 1983, Heggie 1984, Louis & Spurzem 1991). 
Diffusion coefficients are determined self-consistently, including 
the anisotropy of the background distribution of scatterers, but lo- 
cally without any orbit average (Spurzem & Takahashi 1995). Two 
free parameters (A scaling the heat flux, and Xa scaling the colli- 
sional decay of anisotropy) are determined by quantitative compar- 
ison with orbit-averaged Fokker-Planck and direct iV-body mod- 
els, a standard binary heating term due to three-body binary energy 
generation is used; this and the presently used form of the closure 
and all equations can be found in Giersz & Spurzem (1994). Note 
that our present variant (sometimes called a one-flux closure, see 
Louis 1990 and Louis & Spurzem 1991) is the only one, which 
behaves reasonably well for the case of multi-mass systems with 
equipartition and dynamical friction, as was shown in Spurzem & 
Takahashi (1995). We present here an enhancement to the AGM 
method to handle an external tidal field, and also use data of the 
2D FP method worked out by Takahashi (1995, 1996, 1997) and 
iV-body data kindly provided by S. Deiters and D.C. Heggie and E. 
Kim (2003) for comparison purposes. 

The equations of AGM are solved numerically by an im- 
plicit Henyey method as described in Giersz & Spurzem (1994) 
or Spurzem (1996). Collisional terms were evaluated locally with 
self-consistent anisotropic test and background star distributions 
(Giersz & Spurzem 1994, Spurzem & Takahashi 1995). The energy 
generation due to binaries is the same as in the FP models, though 
in the gaseous models it is applied locally, not in an orbit- averaged 
way. 

To describe the process of stars escaping from a cluster in a 
tidal field we adopt simple approximations which are outlined in 
the following. In the spherically symmetric AGM full phase space 
information is reduced to the knowledge of moments of the ve- 
locity distribution up to third order (density, bulk radial mass mo- 
tion, radial and tangential velocity dispersion, radial fluxes of ra- 
dial and tangential energy, see e.g. Spurzem 1994 or Giersz & 
Spurzem 1994 for details of the models). Since the mass and energy 
fluxes drive the quasistatic evolution of the system under relaxation 
timescales, to first order the most significant moments shaping the 
velocity distribution function are density and the radial and tangen- 
tial velocity dispersions. 

In the standard tidal cutoff picture a star is considered to be an 
escaper if its integrals of motion (that of an isolated cluster), energy 
and angular momentum, fulfil certain criteria related to the tidal 
radius r t or tidal cutoff energy E± of the cluster (Takahashi, Lee 
& Inagaki 1997). We approximate the real distribution function of 
the tidally limited cluster in the gaseous model as follows: a tidally 
unlimited local anisotropic Schwarzschild Boltzmann distribution 
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(where p is the mass density, v r ,vg, v$ are the individual stellar ve- 
locities in a local Cartesian coordinate system, whose axes are tan- 
gential to the radial, polar and azimuthal directions on a sphere, u is 
the bulk mass transport velocity, o>, a t are the radial and tangential 
velocity dispersion, respectively) is used to compute the fraction of 
stars X e which would be in the escaper region if a tidal limit would 
be imposed on such a cluster by 
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Here, the index e at the integral denotes an integration over escape 
velocity space; for example in case of an energy criterion we have 
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as escape condition with the potential $ t given at the tidal radius 
(r t = {M/3M g ) 1/3 Rg, where M is the total cluster mass, M G is 
the parent galaxy mass and Rg is the distance between the galaxy 
and cluster centers), v esc is the escape velocity and $(r) is the 
potential at the distance r. For an apocentre criterion, taking into 
account energy and angular momentum conservation (at distance 
r and apocentre distance set to rt), it is straightforward to obtain 
from Eq.0 
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as a condition for escape, where we have used the convenient def- 
initions -u csc ,r = v csc and t4 Clt = v 2 BC • r 2 /(r 2 - r 2 ). The in- 
tegral in Eq. [2] is a 3D-integral over an unbounded domain with 
an ellipsoidal boundary. In order to solve it we first integrate for 
the non-escaping fraction of stars X ne (which are inside the ellip- 
soidal boundary), which is a bound integration domain, and due to 
the normalization of the velocity distribution function we can de- 
termine X e — 1 — X ne . Fortunately, the integral for Xne can be 
solved analytically by using error functions and some similar kind 
of integrals, including Dawson's integral, for which a routine can 
be obtained from Numerical Recipes (Press et al. 1986). In the Ap- 
pendix some more detail of this derivation is given. Here, we show 
the result, using abbreviations: 
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with the definitions of G = a - b and H = b — a . The 
special function T{x) is related to Dawson's integral and defined 
in the Appendix. Similarly to X e we also compute the fraction of 
energy of the stellar system (radial and tangential) belonging to the 
escaper space by 

X r = J v 2 fd 3 v I J v 2 fd 3 v (7) 

e 

Xt = I (vl + vl)fd 3 V j [ (V 2 g + vl)fd 3 V . (8) 
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If X r = Xt = X e all escaping stars would have the same average 
specific energy as the non-escaping ones. Generally this should not 
be the case, since escaping stars would tend to have higher specific 
energies, in such case the difference between X r , Xt and X e tells 
us something about the specific energy of the escapers as compared 
to the non-escaping stars. Again the values are determined by first 
integrating over the bounded domain of the non-escaping stars and 
getting the complement due to the normalization. Our results from 
the Appendix are: 
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where erf, Erf, 1, and J are the error function and specially de- 
fined generalizations of it (see Appendix). 

In a realistic case, however, the distribution function will be 
different from a Schwarzschild-Boltzmann function and the esca- 
per fraction of velocity space will be populated by a few stars only, 
which are on their way out to leave the cluster. Therefore, we as- 
sume that the density of stars p e prone to escape in reality is smaller 
than X e p by a factor k < 1, which will be referred as to filling 
factor. Hence we have p e = kX e p, and for the radial and tan- 
gential energies p re = kX r p r , p te = kXtpt, where p r — pa 2 
and p t = pa 2 . Such procedure can be seen in close connection to 
the ansatz of King's models (King 1966), which just use a lowered 
Maxwellian to model the distribution function of a tidally limited 
cluster. Then our ansatz for mass and energy loss of the cluster is 
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where t CIOSS denotes a crossing time to reach the tidal radius with 
the radial escape velocity, and a is a free parameter with which 
one can describe the unknown process of removal of escaping stars 
from the cluster. 

To complete our model the time evolution of the filling factor 
k has to be described. We are doing this in a close analogy to the 
loss cone description of Frank & Rees (1976) and Amaro-Seoane, 
Freitag & Spurzem (2004) for stars to be swallowed by a central 
black hole. The process which brings stars into the escaper region is 
two-body relaxation, so to the first order we think that the timescale 
ti n to refill the "loss cone" (which is here the escaper region of 
velocity space) is assumed to be 



tin = f3t ry 
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where t rx is the local relaxation time. We keep a free parameter (3 
because some details of the process, e.g. to what extent it is a true 
diffusion process, remain unclear at the moment. The timescale for 
stars to leave the loss-cone is 



Hence we have at each radius r an approximate "diffusion" equa- 
tion describing how stars enter and leave the escaper region of ve- 
locity space: 
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We put the term "diffusion" for this process in quotation marks 
here, because the underlying physical process is angular momen- 
tum diffusion in stellar systems, which is properly described only 
in a 2D Fokker-Planck model, using proper diffusion coefficients. 
The diffusion process tends to establish isotropy, i.e. an equal distri- 
bution of angular momenta - distribution function does not depend 
on angular momentum across the loss-cone and in its vicinity in ve- 
locity space. Here, we model this process in very simplified way by 
the ansatz that the diffusion term is proportional to the difference 
in densities inside (low angular momentum) and outside (high an- 
gular momentum) the loss cone. So, if k = 1, we have an isotropic 
equilibrium and the term vanishes, if k < 1 the diffusion refills the 
loss cone with a rate proportional to 1 /tin, where tin is the standard 
local relaxation time times a factor j3 of order unity. The reader is 
also referred to Amaro-Seoane, Freitag, & Spurzem (2004), where 
a similar concept is used for star accretion onto a massive central 
black hole. 

Suppose, we have to readjust the filling factor k of escaper 
space for some model at time t and radial shell r during the numer- 
ical solution of the gaseous model equations. Then we can consider 
locally p and X e as constant, and find after dividing a factor pX e 
out of the above equation 

dk k 



dt 



to 



+ 
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It can be solved directly only if ti n and t out are held fixed. But we 
use this solution only to advance k(t) from one time step to another, 
then reinitialising it with new values of ij n and t out . If the timestep 
is small enough such that ti n and t out do not change much this is a 
good approximation. 
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Here fco = fc(tci) and Kq = 1 + t ou t/tm; for t — > oo we find 
a stationary solution fcoo = 1/(1 + tin/tout). In our numerical 
models we compute the new filling factor k — k(r) according to 
Eq. ^|in every radial shell between the time steps of the Henyey 
iteration by setting ko = fc(to) and t = to + dt with the model 
time step dt. For small time steps a linear approximation to the 
exponential function would be sufficient, which would again lead 
to linearly discretised form of Eq. lT7l since the computational costs 
are small compared to the Henyey iteration we prefer to be on the 
safe side and use the full exponential function expression Eg. 1181 
for the recomputation of k. 

Note, some special meanings of the parameters a and fi. 
Choosing for example a very small a <C 1 is equivalent to an im- 
mediate removal of escaping stars from the system, as e.g. Chernoff 
& Weinberg (1990) and other FP models usually treat the escapers. 
If a is of the order of one it means that we allow for some time 
before the actual removal of the stars in a complex tidal gravita- 
tional field really takes place. If (3 <C 1 the loss cone region is very 
quickly refilled (practically all the time is full), and if (3 2> 1 the 
loss cone is not refilled. 

This simplified diffusion and escape model, described just by 
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the two time-scales tout and U n , coming with the corresponding 
two parameters a and f3, is being used only in those gaseous model 
shells were E < Et, i.e. the entire shell is still bound to the sys- 
tem and its total energy is smaller than the tidal energy. For shells 
whose energy as a whole is lifted above the tidal energy Et we 
follow a prescription originally proposed by Lee & Ostriker (1987, 
L087). They point out that stars at the tidal energy need very long 
time to actually escape from the cluster, and only if their energy 
is higher than that, they will asymptotically escape with a time 
scale proportional to the crossing time at the tidal radius. The ansatz 
of L087 for the evolution of the phase space distribution function 
/ = f{E, J 2 ) is 

df T /l-E|\ 3 "l 1/2 1 Pitt 

i=- a -n i -(M)J 2^v ~3~ Gpav ■ (i9) 

Note that we have added here the modulus of energies E, E t be- 
cause our sign convention for the energy (E < 0) is different from 
that of L087. p av denotes the average density of the cluster, so the 
last term is actually inversely proportional to a crossing time at the 
tidal radius. Portegies Zwart & Takahashi (1999) find that such a 
model provides a good match between FP models and direct N- 
body simulation. According to their choice of words, there is no 
"crisis" in Fokker-Planck models due to the iV-dependance of the 
dissolution time, provided the parameter ofp is adjusted to a value 
of order unity. L087 varied qfp from 0.2 to 5.0, while Takahashi 
& Portegies Zwart (2000) in an extensive multi-mass study for re- 
alistic globular clusters (with mass spectrum and stellar evolution) 
claim that qfp = 2.5 is the best value. We will discuss the role of 
qfp in our models presented here later; physically qfp brings in 
another time scale of mass loss (rather than the relaxation time), so 
its role is most prominent for small particle numbers, where both 
time scales are comparable. If N is larger, typically 32000 or more, 
the crossing time is small compared to the relaxation time, so the 
latter completely determines the overall evolution of the system. 
But if in the course of the tidal evolution of a star cluster the parti- 
cle number drops, at some moment the terms of L087 will become 
important, and so the time of final dissolution of any cluster will 
depend on qfp- 

In the anisotropic gaseous model we implement the dynamical 
mass loss of L087 simply by applying a mass loss term in the den- 
sities corresponding to Eq |19l if in a radial shell there is E > E t : 

dp r / |£| W /2 1 /4?r 

The tidal energy E t is determined in the gaseous model sim- 
ply by E t = GM(t)/rt(t), where rt results from the standard 
condition that the average density of the entire cluster remains 
constant, and M(t) is the time-dependent total mass (r t (t) = 
(M(t)/Mo) 1/3 n(0), where M and r t (0) are the initial total mass 
and tidal radius, respectively). 

At any time step we redetermine the tidal radius from the 
present mass of the cluster, and remove all shells which fall outside 
of the newly determined tidal radius. This leads typically to a small 
zone where E > E t , but still r < r t in which the L087 dynamical 
mass loss procedure applies. Sometimes, in the very late phases of 
rapid mass loss before dissolution it can happen that shells inside 
the tidal radius have even positive energy (E > 0). If this happens 
the shell is immediately removed, as if it would lie outside the tidal 
radius. 

Our model of mass loss in a tidal field can be determined by in 
total three parameters: a and /3, which describe the time scales of 
the loss cone description for stars entering and leaving the escaper 



Table 1. Parameters of the initial models 



Name 


N 


Uo 


a 


P 


«FP 


k-trh 


AGM-1 


128000 


3 


1 


1 


1 


No - No 


AGM-2 


128000 


6 


1 


1 


1 


No - No 


AGM-3 


128000 


9 


1 


1 


1 


No - No 


AGM-4 


64000 


3 


1 


j 




No - No 


AGM-5 


64000 


6 


1 






No - No 


AGM-6 


64000 


9 


1 


j 


[ 


No - No 


AGM-7 


32000 


3 


1 


1 


1 


No -No 


AGM-8 


32000 


6 


1 


I 


1 


No -No 


AGM-9 


32000 


6 


1 


0.5 


1 


No - Yes 


AGM-10 


32000 


6 


2.5 


0.5 


1 


No - Yes 


AGM-1 1 


32000 


6 


5 


0.5 


1 


No - Yes 


AGM-1 2 


32000 


6 


1 


1 


1 


No - Yes 


AGM-1 3 


32000 


6 


2.5 


1 


1 


No - Yes 


AGM-14 


^ f\f\f\ 

32000 


6 


5 


1 




No - Yes 


AGM-1 5 


32000 


6 


1 


2 


[ 


No - Yes 


AGM-1 6 


32000 


6 


2.5 


2 


1 


No -Yes 


AGM-1 7 


32000 


6 


5 


2 




No - Yes 


AGM-1 8 


32000 


6 


1 


1 





No -No 


AGM-1 9 


32000 


6 


1 


1 


0.2 


No -No 


AGM-20 


32000 


6 


1 


1 


2.5 


No - No 


AGM-2 1 


32000 


6 


1 


1 


5 


No - No 


AGM-22 


32000 


6 


1 


0.5 


1 


Yes - No 


AGM-23 


32000 


6 


1 


1 


1 


Yes - No 


AGM-24 


32000 


6 


1 


2 


1 


Yes - No 


AGM-25 


32000 


9 


1 


1 


1 


No -No 


AGM-26 


16000 


3 


1 


1 


1 


No - No 


AGM-27 


16000 


6 


1 


1 


1 


No - No 


AGM-28 


16000 


6 


1 


2 


1 


No -No 


AGM-29 


16000 


9 


1 


1 


1 


No -No 


AGM-30 


5000 


3 


1 


1 


1 


No - No 


AGM-3 1 


5000 


6 


1 


1 


1 


No - No 


AGM-32 


5000 


6 


1 


1 





No -No 


AGM-33 


5000 


6 


1 


1 


2.5 


No -No 


AGM-34 


5000 


6 


1 


1 


1 


Yes - No 


AGM-35 


5000 


9 


1 


1 


1 


No - No 


AGM-36 


1000 


3 


1 


1 


1 


No - No 


AGM-37 


1000 


6 


1 


1 


1 


No -No 


AGM-38 


1000 


9 


1 


1 


1 


No -No 


IFP-1 


32000 


3 


1 


1 


1 


No -No 


IFP-2 


32000 


6 


1 


1 


1 


No -No 


IFP-3 


1000 


3 


1 


1 


1 


No -No 


IFP-4 


1000 


6 


1 


1 


1 


No -No 


AFP-1 


5000 


6 


1 


1 


1 


No -No 


Nbody-1 


16000 


6 










Nbody-2 


5000 


6 











AGM - anisotropic gaseous models, IFP - isotropic Fokker-Planck models, 
AFP - anisotropic Fokker-Planck models, Nbody - Af-body models. N - 
initial number of stars in the model, Wo - concentration parameter of the 
King model, a and /3 - parameters describing the process of removal of 
escaping stars and replenishing the escape region, respectively (see text - 
discussion of Eqs. 11-16), opp - parameter of Eq. 20 (see text), k - full or 
in equilibrium loss cone of escaping stars (No - equilibrium, Yes - full), trh 
- initial or actual half mass relaxation time used in diffusion equation (No - 
initial, Yes - actual) - see text. 
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loss cone, and Qpp, which scales dynamical mass loss terms at the 
outer boundary. 

Loss terms according to Eqs. I1 1 1 fT2ll 1 3l and l20l are applied as 
additional terms to the gaseous model equations during the Henyey 
iteration. It turns out that both the standard loss terms of Eqs. II II 
I12lll3l as well as the L087 term Eq.[20|are necessary to get a rea- 
sonable fit of mass loss in a cluster in a tidal field as compared 
with direct TV-body as well as FP models. The filling degree k(r) 
can be either set to a constant initial value of k — or k = 1 ac- 
cording to empty or full loss cones, or the stationary value fcoo can 
be used which generally depends on r. For most models discussed 
in this paper we used the stationary value feoo initially, see Table 
1. It is interesting to note that in an unpublished study of Vicary 
(1997) using a Lagrangian gaseous model (Heggie 1984) the same 
qualitative behaviour of mass loss was found as here using a much 
simpler description based upon the loss of mass shells across the 
tidal boundary only. Unfortunately further details of that model are 
not available now. 



3 INITIAL SETUP 

All models discussed in the paper are described by idealized sin- 
gle mass star clusters under the influence of external tidal field. 
The initial positions and velocities of all stars were drawn from a 
King model. The set of initial King models were characterized by 
W = 3, 6 and 9, number of stars by N = 1000, 5000, 16000, 
32000, 64000 and 128000. The range of parameters a, (3 and a FP 
are 1,2.5,5; 0.5,1,2 and 0,0.2,0.5,1,2.5,5, respectively. The 
most natural values for parameters are: a = 1 /3 = 1 and Qpp = 1, 
which means that the time scales for the processes characterized by 
these parameters are described exactly by the local crossing and re- 
laxation times. The initial models are described in Table 1. The data 
for anisotropic Fokker-Planck model (AFP) and A-body model for 
particle number N = 5000 were kindly provided by Kim (2003), 
and TV-body data for N = 16000 by Heggie and Deiters. 

To properly describe, in AGM, the process of mass removal 
connected with loss cone effect and removal of unbounded shells 
it was necessary to increase the spatial resolution of the model, 
particularly close to the tidal boundary. Removal of stars or shells 
from the gaseous model is a big challenge for the numerical algo- 
rithm, when the boundary of the model is closed. Enhanced mass 
loss generates errors of density, velocity, energy and mass distribu- 
tions, which lead to uncontrolled error in the total mass and energy 
of the system. To reduce the error to the acceptable level (a few 
percent) we decided to use, instead of the logarithmically equidis- 
tant grid-points, the mesh, which resolution was enhanced towards 
the tidal boundary. After some numerical experiments the number 
of shells was chosen to: 541, 783, 967 for Wo = 3, 6, 9, respec- 
tively. This guarantee that mass and energy errors was always less 
than 7%. Unfortunately, the code became less efficient. The time of 
calculations is linearly proportional to the number of shells. 

For the computational units the standard TV-body units (Heg- 
gie & Mathieu 1986) were used: total mass M = 1, G = 1 and 
initial total energy of the cluster equal to —1/4. In all presented 
in this paper figures the unit of time is expressed in terms of the 
initial half-mass relaxation time, which for single mass system and 
A-body units is equal to (Spitzer 1987): 



W = 6, N = 32000, o = 1.0 




trh.0 — 



0.138/Vor^ 2 



(21) 



m( 7 /Vo) ' 

where Ao, rno are initial number of stars and half-mass radius 



Figure 1. a: Central density of initial King model Wo = 6 as a function 
of time for gaseous models with 32000 particles, constant a = 1.0 and 
varying j3 = 0.5, 1.0, 2.0 as indicated in the key. Increasing j3 leads mono- 
tonically to increasing lifetimes of the cluster, because the mass loss time 
scale is larger. 



respectively. The value of coefficient in the Coulomb logarithm is 
taken to be 7 = 0.11 (Giersz and Heggie 1994a). 



4 RESULTS 

Fig. la shows the effect of varying /3 on the evolution of the central 
density in the gaseous model with N = 32000 particles in the 
case of constant a = 1. First, one can see that the evolution of the 
system already in pre-collapse varies for different f3\ the reason is 
that the smaller f3 the quicker refilling of the loss cone region and 
the higher mass loss from the system. In other words, if stars are 
removed from the cluster at a certain multiple of the crossing time, 
this process has a different speed relative to relaxation for each N . 
With constant a Fig. lb shows, that a variation of j3 just increases or 
decreases the efficiency of the escape process in a constant amount 
for each particle number. The differences in the dissolution time 
are even greater than in the collapse time. This is connected with 
the fact that the dynamical mass loss described by Eq. 20 is most 
prominent for the small particle numbers - models with smaller j3 
and higher mass loss. 

In contrast to that Figs. 2 show the effect of varying a with 
constant j3. There is almost no effect of changing a on collapse, 
and a small effect on dissolution times. The rate of mass loss is 
only slightly faster for small a than for larger a (see Fig. 2b). This 
is consistent with the picture that the rate of mass loss strongly 
depends on the rate of refilling the loss cone region. Star first have 
to be scattered by the relaxation process to the loss cone region and 
then on the crossing time scale escape from the system. 

From Figs. 1 and 2 one can see that the diffusion time scale 
tin ~ Pt rx , i.e. the parameter has a significant influence on the 
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W n = 6, N = 32000, = 1.0 



W = 6, N = 32000, p = 1.0 




Figure 1 - continued b: As Fig. la, but for the total mass of the system. 



core collapse time and the final dissolution time, by changing the 
mass loss rate. If the mass loss rate is larger, the core collapse accel- 
erates, since the actual half-mass relaxation time becomes smaller 
than the initial one, and obviously the final dissolution time is dif- 
ferent. We have used ctpp = 1.0 for all models, unless stated oth- 
erwise. 

Figs. 3 demonstrate the role of ofp - it hardly changes the 
early evolution, however, the onset of the final dramatic dissolution 
phase can be influenced by it. In the extreme case, where q?fp — > 
the final dissolution is taking extremely long time. This is clearly 
an unphysical case, as has already been stated by L087 and Taka- 
hashi & Portegies Zwart (1998, 2000), Portegies Zwart & Taka- 
hashi (1999). We also show one example, where we have started 
initially with full loss cones, i.e. the phase space part in the escaper 
region is fully populated for all regions of the system. It leads to 
a quick mass loss in the beginning, due to draining of the full loss 
cone in a dynamical time, which cannot be resolved in the figure, 
and thereafter we end up with a somewhat faster evolution. In all 
other models we start with a stationary filling degree of the loss 
cone, koa. 

Now, we turn to comparison of our models with direct iV-body 
results. There are only few models available, and most of them are 
not published, or only partly published. We use here one 16000 N- 
body simulation kindly provided by S. Deiters and D.C. Heggie, for 
Wo = 6, and another one, using 5000 particles, kindly provided by 
E. Kim (2003). 

Comparing these model results with those of our AGM (Figs. 
4) we find that the natural choice of a = 1.0, /3 = 1.0 and 
app = 1.0 provides a fairly good match of the mass loss between 
A^-body and AGM for both particle numbers. However, the core 
collapse times differ by a non-negligible amount. The increase of 
the parameter (5 to values much larger then one will not solve the 
problem, because this will destroy very good agreement for the 
mass loss rate. It should be stressed here the good agreement in 




Figure 2. a: Central density of initial King model Wq = 6 as a function 
of time for gaseous models with 32000 particles, constant f) = 1.0 and 
varying a = 1.0, 2.5, 5.0 as indicated in the key. 



the mass loss rate between A^-body and AGM models. This shows 
the adopted simplified model for the mass loss from the AGM de- 
scribes well the complicated process of mass loss from the stellar 
systems. Finally, it should be stressed that the noisier appearance of 
the Af-body model with 16000 particles as compared to that with 
5000 (see Fig. 4a) is connected with the fact that the larger model 
undergoes gravothermal like oscillations, whereas the smaller one 
shows a smooth post-collapse evolution. However, in the collapse 
phase, as expected, the smaller A^-body model is noisier. 

Now, we show a comparison of our AGM results with 
isotropic ID Fokker-Planck results obtained ky K. Takahashi's 
code. The results are shown in Figs. 5, 6. One sees that now, the 
overall agreement between anisotropic gaseous and isotropic FP 
modes for different King model concentration parameter (Wo = 
3, 6) is quite good for the standard set of parameters: a = 1, 
(3 = 1, ofp = 1, in particular regarding the mass loss, but less 
good for the collapse time. The small adjustment of these parame- 
ters will give a better agreement of the mass loss rate, but still will 
not change much the picture for the collapse time. It seems that gen- 
erally tidally limited anisotropic gaseous models intrinsically show 
a too fast collapse time. AGM predicts a collapse time of about 
60% of the A^-body one. It will be discussed further below in con- 
nection with Fig. 7 that this may be due to an incorrect distribution 
of tidal mass loss across radial shells in our model compared to the 
A-body results (nevertheless we get obvioulsy the total mass loss 
rate well). So, taking for example too much mass out of the core 
would change the collapse time incorrectly. Note that comparing 
individual A^-body models (like we do here) with AGM creates fur- 
ther deviations in core collapse time (see for a single A" = 10000 
A^-body run Spurzem & Aarseth 1996), because, as Giersz & Heg- 
gie (1994a,b) and Giersz & Spurzem (1994) have shown, a proper 
ensemble average of direct A-body models makes the convergence 
to AGM results reliable. In the cited papers we find that AGM pre- 
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W = 6, N = 32000, p = 1.0 W = 6, N = 32000, ct = 1.0, P = 1.0 




5 10 15 20 25 30 35 40 5 10 15 20 25 30 35 40 

"'rhO "'rhO 
Figure 2 - continued b: As Fig. 2a, but for the total mass of the system. Figure 3 - continued b: As Fig. 3a, but for the total mass of the system. 



diets core collapse times in a generally very good agreement with 
iV-body models. 

One might think that a further adjustment of the gaseous 
model parameters (such as the general conductivity A, see e.g. 
Giersz & Spurzem (1994), which scales the collapse time in iso- 
lated models without mass loss) would improve the situation. After 
all, the standard value of A = 0.4977, corresponding to C — 0.104 
in the models of Giersz & Heggie (1994a), has been never checked 
for the case of tidally limited systems. In a larger number of fur- 
ther numerical experiments we find, however, that there is a non- 
trivial relation between the value of A, the core collapse time, and 
the value of /3, which changes the mass loss rate, and thus also the 
core collapse time (see Figs. la,b). It was impossible to find a bet- 
ter combination of A and /3, which both reproduce the mass loss 
and core collapse behaviour more accurately. We concluded that 
the optimal combination of parameters is the use of the standard 
value of A = 0.4977 together with f3 = 1; it provides good agree- 
ment of mass loss rates, but an error in the core collapse times has 
to be tolerated at this point. One possible reason is a still differ- 
ent origin of escapers in the AGM and the TV-body model. Fig. 
shows the cumulative number of escapers as a function of the ra- 
dius at which the last scattering took place lifting them to escape 
energy, comparing AGM and A^-body for 5000 bodies, at a time 
close to the initial model. While both models show that escapers 
can origin from near the core or half-mass radius, it is clear that the 
AGM underestimates the rate of escapers from the halo. The larger 
mass loss from the core the faster the system evolution. We do not 
intend here to further adjust some gaseous models parameters for 
this effect, because this is beyond the scope of this paper. 

One can see from Fig. [SJthat the overall escape rate is very 
well modelled here in the AGM, by looking at the half-mass time 
(time at which the model will contain half of its initial mass) as a 
function of the initial number of stars in comparison to published 
Af-body results. The scaling is proportional to i^ 4 for low N; stan- 



W = 6. N = 32000, a = 1 .0, p = 1 .0 




Figure 3. a: Central density of initial King model Wo = 6 as a function 
of time for gaseous models with 32000 particles, constant /3 = 1.0 and 
a = 1.0. Here we vary the dynamical mass loss constant ctpp as indicated 
in the key. The * symbol denotes a case started with full loss cones, k = 1 
(see text). 
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W„ = 6,a = 1.0,P = 1.0 




W = 6,o = 1.0,P = 1-0 



Figure 4. a: Central density of initial King model Wo = 6 as a function of 
time for gaseous models with 5000 (5k) and 16000 (16k) particles, compar- 
ing data with the corresponding direct Af-body models, see keys. 



dard theory suggests instead a scaling proportional to t r h, if poten- 
tial escapers are immediately removed. Only for large N this scal- 
ing is achieved, while for lower N (below 50.000) it turns over to 
the flatter slope, which is in very good agreement with the theory 
and iV-body results presented by Baumgardt (2001). 

We point out that in the gaseous model gravothermal oscilla- 
tions in a tidally limited, mass-losing system are observed for the 
first times, and they are suppressed by the increasing mass loss to- 
wards the end of the cluster life time. The timestep in the FP model 
was chosen large enough to suppress oscillation for the reasons of 
computational time, so there is no conflict between the results of 
the two methods. Gaseous models are computationally cheap, even 
in the anisotropic tidally limited case. A typical model for Wo = 6 
takes about 1.5 hour CPU time on a Pentium 4 with 3 GHz pro- 
cessor (the CPU time for single mass AGM depends linearly on the 
number of shells, which in turn depends, for tidally limited models, 
only on the initial model concentration, Wo)- So there is no prob- 
lem in resolving the fine structure of the oscillations. Note that in 
cases where there are no post-collapse oscillations (TV = 1000) a 
typical complete run takes only several minutes of CPU time on 
the same computer. In its present form the multi-mass variant of 
AGM would scale in CPU time with r? c , where n c is the number of 
discrete mass components. This is fine for 10 components still, but 
becomes prohibitive at 50 or more components (Giirkan, Freitag & 
Rasio 2004, but note that even in such case we are for large systems 
still much faster than direct A^-body models). Our AGM CPU time 
scales with the third power of n c due to the complete implicit so- 
lution of all equations. It would become much faster if, as usual in 
FP multi-mass studies, the dynamical evolution of each component 
and the interaction via equipartition terms would be decoupled in 
the numerical solution into two steps. Work on this problem is in 
progress. 




Figure 4 - continued b: As Fig. 4a, but for the total mass of the system. 



5 CONCLUSION 

We have derived suitable model equations and their numerical solu- 
tions for star clusters in a tidal field with mass loss in the framework 
of the anisotropic gaseous models (AGM, Louis & Spurzem 1991, 
Spurzem 1994, Giersz & Spurzem 1994). The equations properly 
model the Af-dependence of the mass loss as discussed by Taka- 
hashi & Portegies Zwart (1998, 2000), Portegies Zwart & Taka- 
hashi (1999); since escapers are not removed immediately from the 
system but with a time scale of order of the local crossing time, 
small systems retain their escapers for a longer time with respect to 
their relaxation time scale than large clusters. Our results reproduce 
the evolution of the mass loss of the cluster in good agreement with 
Fokker-Planck (FP) and direct A/-body models; some discrepancy 
in the core collapse time (gaseous model collapses too fast by some 
40%) has to be accepted at the present time. We think that a possi- 
ble reason for this problem is the different radial distribution of the 
origin of escapers in the clusters between AGM and direct A/-body 
model, as seen in Fig. An improvement of this is beyond the 
scope of this paper. Still our gaseous models provide an excellent 
tool to model the global mass loss behaviour with a very efficient 
code (see Fig.[8}. 

Cluster evolution in tidal fields proceeds in three main phases 
according to our results. First, the pre-collapse evolution, which is 
only slightly modified by the mass loss as compared to the isolated 
model. Second, the steady mass loss phase in post-collapse, which 
in the framework of the gaseous model can be described by a sim- 
plified diffusion - escape picture, and drives the system steadily to 
smaller and smaller mass. At some time the mass is small enough 
to lift more and more stars across the tidal energy, and a runaway 
mass loss sets in, which is properly described by the mass loss for- 
mula of Eg. 1201 While in the FP model the diffusion in energy and 
angular momentum is naturally included (which had to be added as 
an additional feature in the gaseous model), the final runaway mass 
loss is not well described by both models, even not in the FP model, 
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Figure 5. a: Central density of initial King model Wo = 3 as a function of 
time for anisotropic gaseous (AGM) and isotropic FP (IFP) models (app = 
1.0) with 1000 (lk) and 32000 (32k) particles, f3 = 1.0 and a = 1.0. 




Figure 5 - continued b: As Fig. 5a, but for the total mass of the system. 



unless it is completed by the term of Eq.^| In the gaseous model it 
is possible to switch off the diffusion across the tidal energy; when 
doing so, we find that the total mass evolution is extremely slow, 
practically halted, even if the mass loss term Eq. 1201 is included. 
Only the presence of steady mass loss caused by diffusion leads to 
conditions where the runaway dissolution of the cluster finally can 
take place. Gaseous models have here again shown their ability as 
excellent tools to analyze the physical processes going on in the 
evolution of heat conducting spheres, as a model for relaxing star 
clusters. 

As shown in Fig. [8] we find that half-mass times (thalf, time 
at which the model will contain half of its initial mass) scale pro- 
portional to £^ 4 . Standard theory predicts to t^alf oc t r h, if poten- 
tial escapers are immediately removed. Our result using the AGM, 
however, is in excellent agreement with an improved theory and 
A^-body results presented by Baumgardt (2001), who points out 
that there is a difference in time between the moment when the 
star reaches an energy necessary to escape and the moment when it 
leaves the cluster. Potential escapers can be scattered back to non- 
escaping energy during this time and become bound again. This 
process is responsible for a deviation from the standard scaling and 

3/4 

leads to than oc t rh for particle numbers at or smaller than 50000. 
This good agreement with A-body results additionally confirms 
the way in which the escapers are treated in AGM and shows that 
AGM can be used with success in simulations of tidally limited star 
cluster. 

For large A^ (32000) we observe gravothermal oscillations 
(Bettwieser & Sugimoto 1984, Goodman 1987, Makino 1996) dur- 
ing the steady mass loss phase in post-collapse. We observe for 
the first time the suppression of post-collapse gravothermal oscilla- 
tions by a critically increasing mass loss at the end of the cluster life 
time (see Fig.|4jt). Note that earlier isotropic FP models of Drukier, 



W = 6,a = 1.0,p = 1.0 




Figure 6. a: Central density of initial King model Wo = 6 as a func- 
tion of time for anisotropic gaseous (AGM) and isotropic FP ( IFP) models 
(a FP = 1.0) with 1000 (lk) and 32000 (32k) particles, (3 = 1.0 and 
a = 1.0. 
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Figure 6 - continued b: As Fig. 6a, but for the total mass of the system. 
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Figure 8. Time at which the total cluster mass is equal to half of the initial 
mass for Wo = 3, Wo = 6 and Wo = 9 as a function of initial number of 
particles AGM models with a = 1, /3 = 1 and a FP = 1.0. 
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r 

Figure 7. Cumulative number of escapers as a function of radius, at which 
the last scattering occurred, which lifted the star to higher than escape en- 
ergy. We compare for 5000 particles, close to the initial Wo = 6 the results 
of AGM and direct A-body. Also the actual core and half-mass radii are 
indicated. 



Fahlman & Richer (1992) and Drukier (1993) also observed cases 
of gravothermal oscillations in tidally limited clusters, but did not 
follow the cluster evolution to full dissolution. The choice of the 
timestep in the FP model shown here was large enough to suppress 
the oscillations. 

The difference between collapse times for AGM and A^-body 
and FP models remains unsolved at the present code version. We 
did extensive parameter studies, many more than actually reported 
in this paper, to find a better optimal combination of parameters 
used to describe the mass loss in AGM (qfp, q, /3), sometimes 
even also varying the standard scaling parameter for the conduc- 
tivity A = 0.4977 (in Heggie's models it is C = 27^k\/IQ, see 
e.g. Giersz & Spurzem (1994). The result is that the standard pa- 
rameters give the best agreement for all other properties between 
AT-body, FP and AGM models, just the core collapse time remains 
too short in AGM. It should be noted that also the agreement be- 
tween FP models and A-body models is not perfect regarding the 
collapse time. We speculate that the radial origin of escapers as de- 
scribed above may be the main reason for the discrepancy. 

Note that our anisotropic gaseous models (AGM) have been 
used in their multi-mass form already with 50 components in 
Giirkan, Freitag, & Rasio (2004) and in Boily et al. (2005) with 15 
components for studies of mass segregation. We want to study the 
future our models with the improved tidal boundary in connection 
with more realistic systems with a stellar mass spectrum (cf. e.g. 
the collaborative experiment, Heggie et al. 1998) and stellar evo- 
lution, and do a further comparison with a set of A'-body models, 
possibly improving also the remaining problems in AGM. 
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APPENDIX 

For clarity we have defined the following integrals in analogy to the 
error function (showed first): 



2 / 2 
erffa:) = — = / exp(— t )dt 
J 
o 

x 

Erf (a:) = -^= / t 2 exp{~t 2 )dt 
J 
o 



2 / 2 

Tlx) — — = / exp(i )dt 



J(x) 



X 

7^1 



t 2 exp(t 2 )dt 



Note the relations 

Erf (a;) = ierf(x) ^= exp(— x 2 ) 

2 y/TY 

1 x 
J(x) = --l(x) + —exp(x 2 ) 

2, yf-K 

2 

T{x) — — = exp(a; )T>(x) 
V 71 " 

where V(x) is Dawson's Integral defined as 



exp(— x ) / t exp(t )dt 



V{x) 



While for the standard error function we use the intrinsic function 
provided by the standard fortran compilers, Dawson's integral has 
to be taken from the Numerical Recipes, Chapter 6.10 (Press et al 
1986). Dawson's Integral vanishes for x — » 0. Since the standard 
method given in the Recipes involves exponential functions this is 
not good for small argument values, therefore for \x\ < 0.2 Daw- 
son's Integral is in the given recipe evaluated by a Taylor series up 
to order x 7 . In our expressions for X e , X r , and X t we have simi- 
larly ill-behaved functions, namely erf (x) /x, Erf (x) / x 3 , T(x) /x, 
and J(x)/x 3 . All these expressions have to be taken for the argu- 
ments G = \Jb 2 — a 2 , b > a or H = \J a? — b 2 , a > b, and they 
approach zero for b — > a. In our numerical computation of these 
functions we therefore also use the following series expansions for 
\x\ < 0.2: 

« ^= f 1 ± -x 2 (l ± —x 2 (l ± —x 2 )) 



X X 

J(x) Eri(x) 
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The ± sign should be taken as a + for I, J (involving oxp(t 2 )) 
and as a — for erf, Erf (involving cxp(— t 2 )). 

While the above Taylor series are completely well behaved 
for x — » 0, it is nevertheless instructive to look at the asymptotic 
forms for X e , X r , and X t which are obtained from the following 
asymptotic equalities for x — > 0: 



T(x) erf (a;) 

x ' x 
J(x) Erf(x) 



exp(±x 2 ) 

V 71 " 

—j=i; exp(±a; 2 ) 
V 71 " 



The use of ± is to be understood as above. With b — a we get the 



results 

2 

X e — erf (a) -= • a • exp(— a ) 

2 2 

X r — erf(a) — • a(l + -a 2 ) ■ exp(— a 2 ) 

Y 7T 3 

X t = erf(a) — • a(l + -a 2 ) • exp(— a 2 ) = X r 

V7T 3 



This same result is obtained approaching b — a from both sides 
(6 < a, b > a), and it reaffirms that in the case of an isotropic 
velocity distribution and equal escape velocities in both the radial 
and tangential direction (i.e. using an energy criterion for escape) 
we have isotropy for the energy of the escaping stars. 
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